---
title: Mechanics
format:
html:
theme: flatly
grid:
body-width: 1000px # Domyślnie jest to ok. 800-900px
margin-width: 250px # Szerokość bocznych paneli (np. TOC)
toc: true
toc-depth: 3
highlight-style: tango
code-line-numbers: true
code-fold: true
code-summary: "Show the code"
code-tools: true
code-block-bg: "rgba(42, 174, 42, 0.02)"
code-block-border-left: "#2aae2a"
code-language-label: true
css: styles.css
math: mathjax
self-contained: true
other-links:
- text: Main page
href: https://dchorazkiewicz.github.io/Mathematics_Physics_Lectures
---
## Description of motion
### Movement in 2D Cartesian coordinates
We could start considering moving particle in 1D but 2D examples are more interesting.
First we have to define the position of the particle in 2D Cartesian coordinates.
The position of the particle is given by the vector $\mathbf{r} = (x, y)$, where $x$ and $y$ are the coordinates of the particle in the $x$ and $y$ axes, respectively. If these coordinates are functions of time, we can write the position vector as $\mathbf{r}(t) = (x(t), y(t))$. This function describes the trajectory of the particle in the $xy$ plane.
$$
[a,b]\rightarrow \mathbb{R}^2:
t\rightarrow \mathbf{r}(t) = (x(t), y(t))
$$
#### Explanation of Symbols:
- $\mathbf{r}$: Position vector of the particle.
- $(x, y)$: Cartesian coordinates representing the particle's position in the $x$ and $y$ axes.
- $\mathbf{r}(t)$: Position vector as a function of time, describing the particle's motion.
- $x(t), y(t)$: Time-dependent functions representing the particle's coordinates in the $x$ and $y$ axes.
- $[a, b]$: Interval of time during which the particle's motion is considered.
- $\mathbb{R}^2$: Two-dimensional Cartesian coordinate space.
- $t$: Time variable, used as a parameter for the position function.
#### How Mathematicians Read This Notation:
Mathematicians interpret the notation as follows:
- "$[a,b] \rightarrow \mathbb{R}^2$" means "a mapping (or function) is defined from the interval $[a, b]$ in time to the two-dimensional Cartesian space $\mathbb{R}^2$."
- "$t \rightarrow \mathbf{r}(t) = (x(t), y(t))$" specifies that for each time $t$ within the interval $[a, b]$, there exists a corresponding position vector $\mathbf{r}(t)$ in $\mathbb{R}^2$, which consists of the time-dependent coordinates $x(t)$ and $y(t)$.
- The overall expression describes a trajectory as a continuous function of time, mapping the progression of time to the corresponding locations in 2D space.
#### Examples
Let us consider two examples of the particle motion in 2D Cartesian coordinates.
$$
\begin{align*}
\mathbf{r}_1(t) &= (t, -t^2 + t) \\
\mathbf{r}_2(t) &= (\cos(t), \sin(t))
\end{align*}
$$
##### Python implementation
```{python}
import numpy as np
import matplotlib.pyplot as plt
def r1(t):
return np.array([t, -t**2 + t])
def r2(t):
return np.array([2 * np.cos(t), 2 * np.sin(t)])
# Define the time ranges
t1 = np.linspace(0, 2, 100)
t2 = np.linspace(0, 2 * np.pi, 100)
# Create a side-by-side plot layout
fig, axes = plt.subplots(nrows=1, ncols=2, figsize=(12, 5)) # One row, two columns
# First plot
axes[0].plot(r1(t1)[0], r1(t1)[1])
axes[0].set_xlabel('x')
axes[0].set_ylabel('y')
axes[0].set_title('Plot 1: r1(t)')
# Second plot
axes[1].plot(r2(t2)[0], r2(t2)[1])
axes[1].set_xlabel('x')
axes[1].set_ylabel('y')
axes[1].set_title('Plot 2: r2(t)')
# Adjust aspect ratios if needed
axes[0].set_aspect('equal', 'box')
axes[1].set_aspect('equal', 'box')
# Optimize layout
plt.tight_layout()
plt.show()
```
##### HTML implementation
<iframe src="html_anim/mechanics/kinematics_1_position_vector.html" width="100%" height="900px" frameborder="0"></iframe>
[html_anim/mechanics/kinematics_1_position_vector](html_anim/mechanics/kinematics_1_position_vector.html)
## Velocity
The velocity of the particle is the derivative of the position vector with respect to time. The velocity vector is defined as
$$
\mathbf{v}(t) = \frac{d\mathbf{r}(t)}{dt} = \left(\frac{dx(t)}{dt}, \frac{dy(t)}{dt}\right)
$$
The velocity vector describes the speed and direction of the particle at any time $t$. The magnitude of the velocity vector is the speed of the particle. The direction of the velocity vector is the direction of motion of the particle.
```{python}
import numpy as np
import matplotlib.pyplot as plt
def r(t):
return np.array([t, -t**2 + t])
def v(t):
return np.array([1, -2*t + 1])
t = np.linspace(0, 10, 100)
fig, axes = plt.subplots(figsize=(8, 6))
# whole curve
axes.plot(r(t)[0], r(t)[1], label='r(t)')
# velocity vector at t=3
axes.quiver(r(3)[0], r(3)[1], v(3)[0], v(3)[1], angles='xy', scale_units='xy', scale=1, color='red', label='v(3)')
axes.text(r(3)[0] + v(3)[0], r(3)[1] + v(3)[1], 'v(3)', color='red')
#velocity vector at t=6
axes.quiver(r(6)[0], r(6)[1], v(6)[0], v(6)[1], angles='xy', scale_units='xy', scale=1, color='blue', label='v(6)')
axes.text(r(6)[0] + v(6)[0], r(6)[1] + v(6)[1], 'v(6)', color='blue')
axes.set_xlabel('x')
axes.set_ylabel('y')
axes.legend()
plt.show()
```
<iframe src="html_anim/mechanics/kinematics_1_velocity.html" width="100%" height="900px" frameborder="0"></iframe>
[html_anim/mechanics/kinematics_1_velocity](html_anim/mechanics/kinematics_1_velocity.html)
## Acceleration
The acceleration of the particle is the derivative of the velocity vector with respect to time. The acceleration vector is defined as
$$
\mathbf{a}(t) =
\frac{d\mathbf{v}(t)}{dt} =
\left(\frac{dv_x(t)}{dt}, \frac{dv_y(t)}{dt}\right) =
\left(\frac{d^2x(t)}{dt^2}, \frac{d^2y(t)}{dt^2}\right)
$$
The acceleration vector describes the rate of change of the velocity of the particle at any time $t$. The magnitude of the acceleration vector is the rate of change of the speed of the particle. The direction of the acceleration vector is the direction of the change of the velocity of the particle.
```{python}
import numpy as np
import matplotlib.pyplot as plt
def r(t):
return np.array([t, -t**2 + t])
def v(t):
return np.array([1, -2*t + 1])
def a(t):
return np.array([0, -2])
t = np.linspace(0, 4, 100)
fig, axes = plt.subplots(figsize=(8, 6))
# whole curve
axes.plot(r(t)[0], r(t)[1], label='r(t)')
# velocity vector at t=3
axes.quiver(r(3)[0], r(3)[1], v(3)[0], v(3)[1], angles='xy', scale_units='xy', scale=1, color='red', label='v(3)')
axes.text(r(3)[0] + v(3)[0], r(3)[1] + v(3)[1], 'v(3)', color='red')
# acceleration vector at t=3
axes.quiver(r(3)[0], r(3)[1], a(3)[0], a(3)[1], angles='xy', scale_units='xy', scale=1, color='blue', label='a(3)')
axes.text(r(3)[0] + a(3)[0], r(3)[1] + a(3)[1], 'a(3)', color='blue')
axes.set_xlabel('x')
axes.set_ylabel('y')
axes.legend()
plt.show()
```
<iframe src="html_anim/mechanics/kinematics_1_acc.html" width="100%" height="900px" frameborder="0"></iframe>
[html_anim/mechanics/kinematics_1_acc](html_anim/mechanics/kinematics_1_acc.html)
### Numerical approximation of velocity and acceleration
<iframe src="html_anim/mechanics/kinematics_numerical.html" width="100%" height="900px" frameborder="0"></iframe>
[html_anim/mechanics/kinematics_numerical](html_anim/mechanics/kinematics_numerical.html)
## Newton's laws of motion
Revolutions in physics started with Newton's laws of motion. These laws describe the motion of particles in space. The first law states that a particle moves with constant velocity if no external forces act on it. The second law states that the acceleration of a particle is proportional to the force acting on it. The third law states that forces always occur in pairs. If one object exerts a force on another object, the second object exerts an equal and opposite force on the first object. These laws are the foundation of classical mechanics.
### Newton's first law
Newton's first law states that a particle moves with constant velocity if no external forces act on it. This law is also known as the law of inertia. The law of inertia states that an object at rest stays at rest and an object in motion stays in motion with the same speed and in the same direction unless acted upon by an external force. This law is a consequence of the conservation of momentum. The momentum of a particle is the product of its mass and velocity. The momentum of a particle is conserved if no external forces act on it.
### Newton's second law
Newton's second law states that the acceleration of a particle is proportional to the force acting on it. The acceleration of a particle is the rate of change of its velocity. The force acting on a particle is the product of its mass and acceleration. The force acting on a particle is equal to the rate of change of its momentum. The second law can be written as
$$
\mathbf{F}(x, y, z, t) = m \mathbf{a}(x, y, z, t)
$$
where $\mathbf{F}(t)=(F_x(t), F_y(t), F_z(t))$ is the force acting on the particle, $m$ is the mass of the particle, and
$\mathbf{a}=(x''(t), y''(t),z''(t))$ is the acceleration of the particle.
Above equation can be written as a set of three equations
$$
\begin{align*}
\frac{d^2 x(t)}{dt^2} &= \frac{F_x(x, y, z, t)}{m} \\
\frac{d^2 y(t)}{dt^2} &= \frac{F_y(x, y, z, t)}{m}\\
\frac{d^2 z(t)}{dt^2} &= \frac{F_z(x, y, z, t)}{m}
\end{align*}
$$
#### Bounduary conditions
To solve these equations we need to know the initial position and velocity of the particle. These are the boundary conditions of the problem. The boundary conditions are the initial values of the position and velocity of the particle. The boundary conditions determine the trajectory of the particle in space. The boundary conditions can be used to solve the differential equations of motion.
### Newton's third law
Newton's third law states that forces always occur in pairs. If one object exerts a force on another object, the second object exerts an equal and opposite force on the first object. This law is also known as the law of action and reaction. The law of action and reaction states that for every action there is an equal and opposite reaction. This law is a consequence of the conservation of momentum. The momentum of a system is conserved if no external forces act on it.
### Solving the Second-Order Equation Iteratively
#### 1. Rewriting the Second-Order Equation
The second-order equation is:
$$
\frac{d^2 x(t)}{dt^2} = \frac{F_x(x, y, z, t)}{m}.
$$
We introduce $v_x(t)$, the velocity in the $x$-direction, such that:
$$
v_x(t) = \frac{dx(t)}{dt}.
$$
This converts the equation into a system of two first-order equations:
$$
\begin{align}
\frac{dx(t)}{dt} &= v_x(t), \\
&\\
\frac{dv_x(t)}{dt} &= \frac{F_x(x, y, z, t)}{m}.
\end{align}
$$
#### 2. Iterative Procedure
We solve these equations step-by-step using numerical methods. For simplicity, let's use **Euler's method**:
- Let $x_n = x(t_n)$ and $v_{x,n} = v_x(t_n)$.
- Given a small time step $\Delta t$, the update rules are:
$$
\begin{align}
x_{n+1} &= x_n + v_{x,n} \Delta t, \\
v_{x,n+1} &= v_{x,n} + \frac{F_x(x_n, y_n, z_n, t_n)}{m} \Delta t.
\end{align}
$$
We iterate these equations to compute the trajectory.
#### 3. Focus on the $x$-Component
To make this concrete, let's assume:
- A force $F_x(x, y, z, t) = -kx^3$, where $k$ is a constant.
- We ignore $y$ and $z$ for simplicity (focus on $x$-component only).
This gives:
$$
\frac{dv_x(t)}{dt} = \frac{-k (x(t))^3}{m}.
$$
The iterative update rules become:
$$
\begin{cases}
x_{n+1} = x_n + v_{x,n} \Delta t, \\
v_{x,n+1} = v_{x,n} - \frac{k x_n^3}{m} \Delta t.
\end{cases}
$$
#### 4. Numerical Implementation
Here's how we can implement this:
```{python}
import numpy as np
import matplotlib.pyplot as plt
# Parameters
k = 1.0 # Spring constant
m = 1.0 # Mass
dt = 0.01 # Time step
steps = 1000 # Number of steps
x0 = 1.0 # Initial position
v0 = 0.0 # Initial velocity
def force(x):
return -k * x*x*x
# Arrays to store time, position, and velocity
time = np.linspace(0, steps * dt, steps)
x = np.zeros(steps)
v = np.zeros(steps)
# Initial conditions
x[0] = x0
v[0] = v0
# Iterative solution using Euler's method
for n in range(steps - 1):
# Update position
x[n + 1] = x[n] + v[n] * dt
# Update velocity
v[n + 1] = v[n] + force(x[n]) / m * dt
# Plot results
plt.figure(figsize=(10, 5))
plt.plot(time, x, label='Position (x)')
plt.plot(time, v, label='Velocity (v)')
plt.xlabel('Time (t)')
plt.ylabel('Position/Velocity')
plt.title('Harmonic Oscillator Solution')
plt.legend()
plt.grid()
plt.show()
```
## Basic examples
### Projectile motion
Equations of motion for a particle in free fall are
$$
\begin{align*}
\frac{d^2 x(t)}{dt^2} &= 0 \\
\frac{d^2 y(t)}{dt^2} &= -g
\end{align*}
$$
where $g$ is the acceleration due to gravity. The acceleration of the particle in the $x$ direction is zero. The acceleration of the particle in the $y$ direction is equal to the acceleration due to gravity. The force acting on the particle is the force of gravity.
#### Analytical solution
To find the position of the particle as a function of time, we integrate the equations of motion with respect to time and apply the initial conditions.
**1. Horizontal motion ($x$-direction)**
First, we integrate the acceleration $\frac{d^2 x}{dt^2} = 0$ to find the velocity:
$$
v_x(t) = \int 0 \, dt = C_1
$$
At time $t=0$, the initial velocity is $v_{0x}$. Substituting this into the equation gives $C_1 = v_{0x}$. Thus:
$$
v_x(t) = v_{0x}
$$
Next, we integrate the velocity to find the position:
$$
x(t) = \int v_{0x} \, dt = v_{0x}t + C_2
$$
At time $t=0$, the initial position is $x_0$. Substituting this gives $C_2 = x_0$. The final equation for horizontal position is:
$$
x(t) = v_{0x} t + x_0
$$
**2. Vertical motion ($y$-direction)**
First, we integrate the acceleration $\frac{d^2 y}{dt^2} = -g$ to find the velocity:
$$
v_y(t) = \int -g \, dt = -gt + C_3
$$
At time $t=0$, the initial velocity is $v_{0y}$. Substituting this gives $C_3 = v_{0y}$. Thus:
$$
v_y(t) = -gt + v_{0y}
$$
Next, we integrate the velocity to find the position:
$$
y(t) = \int (-gt + v_{0y}) \, dt = -\frac{1}{2}gt^2 + v_{0y}t + C_4
$$
At time $t=0$, the initial position is $y_0$. Substituting this gives $C_4 = y_0$. The final equation for vertical position is:
$$
y(t) = -\frac{1}{2} g t^2 + v_{0y} t + y_0
$$
**Summary**
The complete analytical solution describing the trajectory of the particle is:
$$
\begin{align*}
x(t) &= v_{0x} t + x_0 \\
y(t) &= -\frac{1}{2} g t^2 + v_{0y} t + y_0
\end{align*}
$$
#### Numerical solution
```{python}
import numpy as np
import matplotlib.pyplot as plt
# Constants
F = (0, -9.81) # Force of gravity (N)
m = 1 # Mass of the particle (kg)
x_0, y_0 = 0, 100 # Initial position (m)
v_x0, v_y0 = 1, 0 # Initial velocity (m/s)
def simulate_trajectory(F, m, x_0, y_0, v_x0, v_y0, t_max, steps):
# Time discretization
t = np.linspace(0, t_max, steps)
h = t[1] - t[0] # Time step
# Initialize position and velocity
x = [x_0]
y = [y_0]
v_x, v_y = v_x0, v_y0
# Euler integration
for i in range(1, len(t)):
a_x, a_y = F[0] / m, F[1] / m # Acceleration
v_x += a_x * h # Update velocity
v_y += a_y * h
x_next = x[-1] + v_x * h # Update position
y_next = y[-1] + v_y * h
# Stop if the particle hits the ground
if y_next < 0:
break
x.append(x_next)
y.append(y_next)
return x, y, t[:len(x)] # Return trajectory and corresponding time
# Simulate the trajectory
x, y, t = simulate_trajectory(F, m, x_0, y_0, v_x0, v_y0, t_max=5, steps=100)
# Create grid for vector field
X, Y = np.meshgrid(np.linspace(-1, 6, 10), np.linspace(0, 120, 5))
U = np.zeros_like(X) # Horizontal component of g is 0
V = -9.81 * np.ones_like(Y) # Vertical component of g is constant
# Time to display the ball and force
t_display = 2.0 # Time at which to show the ball and force
idx = np.argmin(np.abs(t - t_display)) # Find the closest index for the given time
# Visualization
fig, ax = plt.subplots(figsize=(10, 6))
ax.plot(x, y, label="Trajectory")
ax.axhline(0, color='green', linestyle='dashed', label='Ground level')
ax.scatter(x_0, y_0, color='red', label='Initial position')
# Add vector field
ax.quiver(X, Y, U, V, color='blue', alpha=0.3, scale=200
, width=0.002,label='Gravitational field')
# Add the ball at the selected time
ball_x, ball_y = x[idx], y[idx]
ax.scatter(ball_x, ball_y, color='orange', s=100, label='Ball (t=2s)')
# Add the force vector acting on the ball
force_x, force_y = F[0] / m, F[1] / m # Force per unit mass
ax.quiver(ball_x, ball_y, force_x, force_y, color='red', angles='xy', scale_units='xy', scale=0.5, label='Force on Ball')
# Labels and title
ax.set_xlabel('x (m)')
ax.set_ylabel('y (m)')
ax.set_title('Particle Trajectory with Gravitational Field and Force on Ball')
ax.legend()
plt.show()
```
<iframe src="html_anim/mechanics/dynamics_projectile_motion.html" width="100%" height="700px" frameborder="0"></iframe>
[html_anim/mechanics/dynamics_projectile_motion](html_anim/mechanics/dynamics_projectile_motion.html)
### Harmonic oscillation
Let $\mathbb{F}=-k x$ be the force acting on the particle. The equation of motion is
$$
\frac{d^2 x(t)}{dt^2} = -\frac{k}{m} x(t)
$$
where $k$ is the spring constant and $m$ is the mass of the particle. The acceleration of the particle is proportional to its position, and the force acting on the particle is the force of a spring, equal to the spring constant times the position of the particle.
```{python}
import numpy as np
import matplotlib.pyplot as plt
# Updated constant
k = 0.5 # Updated spring constant
displacements = [1, 2, 3, 4, 5] # Displacements (x)
forces = [-k * x for x in displacements] # Corresponding forces (F)
# Positions for visualization
y_positions = np.linspace(1, len(displacements) + 1, len(displacements)) # Avoid overlapping
# Visualization
fig, ax = plt.subplots(figsize=(8, 6))
# Draw force vectors for each displacement
for x, F, y in zip(displacements, forces, y_positions):
# Draw ball position
ax.scatter(x, y, color='orange', s=100, label=f"x = {x}, F = {F:.2f}" if y == y_positions[0] else "")
# Draw force vector
ax.quiver(x, y, F, 0, angles='xy', scale_units='xy', scale=1, color='blue', width=0.005)
# Add annotation for calculations
ax.text(x + F / 2, y + 0.2, f"F = {-k:.1f}*{x} = {F:.2f}", fontsize=9, color='black', alpha=.8)
# Labels and title
ax.axhline(0, color='black', linewidth=0.5, linestyle='dashed') # Equilibrium line
ax.set_xlim(0, 6)
ax.set_ylim(0, len(displacements) + 2)
ax.set_xlabel("Displacement $x$")
ax.set_ylabel("Vertical position (for clarity)")
ax.set_title("Force Linearly Dependent on Displacement with Annotations")
plt.grid(True, linestyle='--', alpha=0.6)
plt.show()
```
#### Analytical solution
The analytical solution of the equation of motion for a particle in simple harmonic motion is
$$
x(t) = A \sin(\omega t) + B \cos(\omega t)
$$
where $A$ and $B$ are the amplitudes of the particle, and $\omega$ is the angular frequency of the particle
$$
\omega = \sqrt{\frac{k}{m}}
$$
The angular frequency is equal to the square root of the spring constant divided by the mass of the particle. The phase angle $\phi$ determines the initial phase of the particle.
#### Numerical solution
```{python}
import numpy as np
import matplotlib.pyplot as plt
# Parameters
k = 1 # Spring constant
m = 1 # Mass of the particle
x_0 = 1 # Initial position
v_0 = 0 # Initial velocity
# Derived parameters
omega = np.sqrt(k / m) # Angular frequency
# Time array
t = np.linspace(0, 50, 500)
dt = t[1] - t[0]
# Numerical solution (Euler method)
x_num = [x_0]
v = v_0
for i in range(1, len(t)):
a = -k / m * x_num[-1] # Acceleration
v += a * dt # Update velocity
x_num.append(x_num[-1] + v * dt) # Update position
# Analytical solution
A = x_0
B = v_0 / omega
x_analytical = B * np.sin(omega * t) + A * np.cos(omega * t)
# Plotting the position vs time
plt.figure(figsize=(8, 4))
plt.plot(t, x_num, label="Numerical Solution", linestyle='--')
plt.plot(t, x_analytical, label="Analytical Solution", linestyle='-')
plt.xlabel('Time (t)')
plt.ylabel('Position (x)')
plt.title('Harmonic Motion in One Dimension')
plt.legend()
plt.grid()
plt.show()
```
<iframe src="html_anim/mechanics/dynamics_harmonic_oscilator.html" width="100%" height="700px" frameborder="0"></iframe>
[html_anim/mechanics/dynamics_harmonic_oscilator](html_anim/mechanics/dynamics_harmonic_oscilator.html)
### Pendulum Motion
#### Introduction
The motion of a simple pendulum is a classic example of harmonic motion. A pendulum consists of a mass $m$ (called the bob) attached to a string or rod of length $L$, which is fixed at one end and free to swing back and forth under the influence of gravity. For small angles, the motion can be approximated as simple harmonic motion.
#### Equations of Motion
The equation of motion for a pendulum is derived from Newton's second law. The force acting on the pendulum is the component of the gravitational force tangential to the arc of its swing:
$$
F = -mg \sin(\theta)
$$
```{python}
import numpy as np
import matplotlib.pyplot as plt
# Constants
pivot = (0, 10) # Point of attachment of the pendulum
L = 10 # Length of the pendulum
theta1 = np.radians(30) # First angle of displacement (in radians)
theta2 = np.radians(60) # Second angle of displacement (in radians)
g = 9.81 # Gravitational acceleration
m = 0.5 # Mass of the pendulum bob
# Calculate bob position for both angles
x_bob1 = pivot[0] + L * np.sin(theta1) # x-coordinate for θ=30°
y_bob1 = pivot[1] - L * np.cos(theta1) # y-coordinate for θ=30°
x_bob2 = pivot[0] + L * np.sin(theta2) # x-coordinate for θ=60°
y_bob2 = pivot[1] - L * np.cos(theta2) # y-coordinate for θ=60°
# Calculate forces for both angles
F_gravity = m * g # Magnitude of gravitational force
F_tangential1 = F_gravity * np.sin(theta1) # Tangential force for θ=30°
F_tangential2 = F_gravity * np.sin(theta2) # Tangential force for θ=60°
# Circle parameters for light shading
circle = plt.Circle(pivot, L, color='gray', fill=False, linestyle='dashed', alpha=0.5)
# Visualization setup
fig, ax = plt.subplots(figsize=(8, 8))
# Add circular path for the pendulum (single circle for both cases)
ax.add_artist(circle) # Add the circular path
# Draw pendulum lines for both cases
ax.plot([pivot[0], x_bob1], [pivot[1], y_bob1], color='black', linewidth=1.5, label='Pendulum (θ=30°)')
ax.plot([pivot[0], x_bob2], [pivot[1], y_bob2], color='black', linewidth=1.5, linestyle='dotted', label='Pendulum (θ=60°)')
# Draw bobs as larger circles for both cases
bob1 = plt.Circle((x_bob1, y_bob1), 0.5, color='orange', zorder=5) # Bob for θ=30°
bob2 = plt.Circle((x_bob2, y_bob2), 0.5, color='green', zorder=5) # Bob for θ=60°
ax.add_artist(bob1)
ax.add_artist(bob2)
# Draw force vectors for θ=30°
ax.quiver(x_bob1, y_bob1, 0, -F_gravity, angles='xy', scale_units='xy', scale=1, color='red', label='Gravitational Force (mg)')
ax.quiver(
x_bob1, y_bob1,
-F_tangential1 * np.cos(theta1), -F_tangential1 * np.sin(theta1),
angles='xy', scale_units='xy', scale=1, color='blue', label=r'Tangential Force ($-mg \sin(\theta)$)'
)
# Draw force vectors for θ=60°
ax.quiver(x_bob2, y_bob2, 0, -F_gravity, angles='xy', scale_units='xy', scale=1, color='darkred')
ax.quiver(
x_bob2, y_bob2,
-F_tangential2 * np.cos(theta2), -F_tangential2 * np.sin(theta2),
angles='xy', scale_units='xy', scale=1, color='darkblue', label=r'Tangential Force ($-mg \sin(\theta)$, θ=60°)'
)
# Labels and title
ax.set_xlim(-L - 2, L + 2)
ax.set_ylim(-5, 12)
ax.set_aspect('equal', adjustable='box')
ax.axhline(pivot[1], color='black', linestyle='dashed', linewidth=1.5) # Horizontal reference
ax.set_title("Forces Acting on a Pendulum (Two Positions)")
ax.set_xlabel("Horizontal Position")
ax.set_ylabel("Vertical Position")
ax.legend()
plt.grid(True, linestyle='--', alpha=0.6)
plt.show()
```
Using the relationship $F = ma$ and the angular acceleration $a = L \frac{d^2 \theta}{dt^2}$, we get:
$$
L \frac{d^2 \theta}{dt^2} = -g \sin(\theta)
$$
Dividing through by $L$, the equation becomes:
$$
\frac{d^2 \theta}{dt^2} + \frac{g}{L} \sin(\theta) = 0
$$
For small angles ($\sin(\theta) \approx \theta$), the equation simplifies to:
$$
\frac{d^2 \theta}{dt^2} + \frac{g}{L} \theta = 0
$$
This is the equation for simple harmonic motion with angular frequency:
$$
\omega = \sqrt{\frac{g}{L}}
$$
<iframe src="html_anim/mechanics/dynamics_pendulum.html" width="100%" height="700px" frameborder="0"></iframe>
[html_anim/mechanics/dynamics_pendulum](html_anim/mechanics/dynamics_pendulum.html)
#### Analytical Solution
The analytical solution for small angles is:
$$
\theta(t) = \theta_0 \cos(\omega t + \phi)
$$
where:
- $\theta_0$ is the initial angular displacement.
- $\phi$ is the phase constant, determined by initial conditions.
#### Numerical Solution
For larger angles, the small-angle approximation does not hold, and we need to solve the original nonlinear equation numerically. We can use the finite difference method to approximate the solution.
```{python}
import numpy as np
import matplotlib.pyplot as plt
# Parameters
g = 9.81 # Acceleration due to gravity (m/s^2)
L = 1.0 # Length of the pendulum (m)
theta_0 = np.pi / 4 # Initial angle (radians)
# Time settings
T = 10 # Total time (s)
N = 10_000 # Number of time steps
dt = T / N
# Arrays to store values
t = np.linspace(0, T, N)
theta_full = np.zeros(N)
omega_full = np.zeros(N) # Angular velocity for full sin(theta)
theta_approx = np.zeros(N)
omega_approx = np.zeros(N) # Angular velocity for sin(theta) ~ theta
# Initial conditions
theta_full[0] = theta_0
omega_full[0] = 0
theta_approx[0] = theta_0
omega_approx[0] = 0
# Numerical integration (Euler method)
for i in range(1, N):
# Full sin(theta)
alpha_full = -(g / L) * np.sin(theta_full[i - 1])
omega_full[i] = omega_full[i - 1] + alpha_full * dt
theta_full[i] = theta_full[i - 1] + omega_full[i] * dt
# Approximation sin(theta) ~ theta
alpha_approx = -(g / L) * theta_approx[i - 1]
omega_approx[i] = omega_approx[i - 1] + alpha_approx * dt
theta_approx[i] = theta_approx[i - 1] + omega_approx[i] * dt
# Plotting
plt.figure(figsize=(12, 6))
plt.plot(t, theta_full, label="Full sin(theta)", color="blue")
plt.plot(t, theta_approx, label="sin(theta) ~ theta", color="red", linestyle="--")
plt.title("Pendulum Motion: Full sin(theta) vs Approximation")
plt.xlabel("Time (s)")
plt.ylabel("Angle (radians)")
plt.grid()
plt.legend()
plt.show()
```
<iframe src="html_anim/mechanics/dynamics_pendulum_math_real_numerical.html" width="100%" height="700px" frameborder="0"></iframe>
[html_anim/mechanics/dynamics_pendulum_math_real_numerical](html_anim/mechanics/dynamics_pendulum_math_real_numerical.html)
#### Energy Analysis
The total mechanical energy of the pendulum is conserved (assuming no damping). The total energy is the sum of kinetic and potential energy:
$$
E = K + U
$$
- Kinetic energy: $K = \frac{1}{2} m v^2$
- Potential energy: $U = m g h$
For a pendulum:
$$
v = L \frac{d\theta}{dt}, \quad h = L (1 - \cos\theta)
$$
Thus, the total energy becomes:
$$
E = \frac{1}{2} m L^2 \left(\frac{d\theta}{dt}\right)^2 + m g L (1 - \cos\theta)
$$
The energy remains constant over time, which can be verified numerically.
```{python}
# Calculate energies
m=1
kinetic_energy = 0.5 * m * (L * omega_full)**2
potential_energy = m * g * L * (1 - np.cos(theta_full))
total_energy = kinetic_energy + potential_energy
# Plot energies
plt.figure(figsize=(10, 6))
plt.plot(t, kinetic_energy, label="Kinetic Energy")
plt.plot(t, potential_energy, label="Potential Energy")
plt.plot(t, total_energy, label="Total Energy", linestyle="dashed")
plt.title("Energy of the Pendulum")
plt.xlabel("Time (s)")
plt.ylabel("Energy (J)")
plt.grid()
plt.legend()
plt.show()
```
### Circular Motion
#### Introduction
Circular motion refers to the movement of an object along a circular path. This motion can be uniform (constant speed) or non-uniform (variable speed). For simplicity, we will consider uniform circular motion, where the speed of the object is constant. The position, velocity, and acceleration vectors in circular motion exhibit interesting properties:
- The velocity vector is always tangent to the circle.
- The acceleration vector (centripetal acceleration) always points toward the center of the circle.
#### Equations of Motion
Assume an object moves along a circular path of radius $R$ with constant angular velocity $\omega$.
##### Position Vector
The position of the object at time $t$ can be described in terms of the angle $\theta(t)$ it makes with the reference axis:
$$
\mathbf{r}(t) = R \cos(\omega t) \hat{i} + R \sin(\omega t) \hat{j}
$$
#### Velocity Vector
The velocity is the derivative of the position vector with respect to time:
$$
\mathbf{v}(t) = \frac{d\mathbf{r}(t)}{dt} = -R \omega \sin(\omega t) \hat{i} + R \omega \cos(\omega t) \hat{j}
$$
The magnitude of the velocity is constant and equals:
$$
|\mathbf{v}(t)| = R \omega
$$
**Velocity and Position Perpendicularity:**
The dot product of the position vector $\mathbf{r}(t)$ and velocity vector $\mathbf{v}(t)$ is:
$$
\mathbf{r}(t) \cdot \mathbf{v}(t) = \left[R \cos(\omega t)\right] \left[-R \omega \sin(\omega t)\right] + \left[R \sin(\omega t)\right] \left[R \omega \cos(\omega t)\right]
$$
Simplifying:
$$
\mathbf{r}(t) \cdot \mathbf{v}(t) = -R^2 \omega \cos(\omega t) \sin(\omega t) + R^2 \omega \sin(\omega t) \cos(\omega t) = 0
$$
This confirms that $\mathbf{r}(t)$ and $\mathbf{v}(t)$ are perpendicular.
#### Acceleration Vector
The acceleration is the derivative of the velocity vector with respect to time:
$$
\mathbf{a}(t) = \frac{d\mathbf{v}(t)}{dt} = -R \omega^2 \cos(\omega t) \hat{i} - R \omega^2 \sin(\omega t) \hat{j} = -R \omega^2 \mathbf{r}(t)
$$
The acceleration vector always points toward the center of the circle (centripetal acceleration), and its magnitude is:
$$
|\mathbf{a}(t)| = R \omega^2 = \frac{v^2}{R}
$$
Last relation is very important because it shows that the acceleration is proportional to the square of the velocity and inversely proportional to the radius of the circle.
This result will we very useful in the next section when we will discuss the gravity.
#### Numerical Example
We will visualize the position, velocity, and acceleration at two distinct points along the circular path.
```{python}
import numpy as np
import matplotlib.pyplot as plt
# Parameters
R = 1 # Radius of the circle
omega = np.pi/4 # Angular velocity (rad/s)
T = 1 # Time period (s)
t_points = [0, T] # Time points to analyze
# Calculate vectors at specific points
positions = [(R * np.cos(omega * t), R * np.sin(omega * t)) for t in t_points]
velocities = [(-R * omega * np.sin(omega * t), R * omega * np.cos(omega * t)) for t in t_points]
accelerations = [(-R * omega**2 * np.cos(omega * t), -R * omega**2 * np.sin(omega * t)) for t in t_points]
# Plot the circle and vectors
fig, ax = plt.subplots(figsize=(8, 8))
theta = np.linspace(0, 2 * np.pi, 100)
x_circle = R * np.cos(theta)
y_circle = R * np.sin(theta)
ax.plot(x_circle, y_circle, label="Circular Path")
for i, t in enumerate(t_points):
x, y = positions[i]
vx, vy = velocities[i]
axx, axy = accelerations[i]
ax.quiver(x, y, vx, vy, color="red", angles="xy", scale_units="xy", scale=1, label=f"Velocity at t={t:.2f}s" if i == 0 else None)
ax.quiver(x, y, axx, axy, color="blue", angles="xy", scale_units="xy", scale=1, label=f"Acceleration at t={t:.2f}s" if i == 0 else None)
ax.set_aspect('equal', adjustable='box')
ax.set_title("Circular Motion: Velocity and Acceleration")
ax.set_xlabel("x")
ax.set_ylabel("y")
ax.set_xlim(-2, 2)
ax.set_ylim(-2, 2)
ax.legend()
ax.grid()
plt.show()
```
## Gravity
### Introduction
The Universal Law of Gravitation, formulated by Sir Isaac Newton, describes the gravitational force between two masses. It states that every particle of matter in the universe attracts every other particle with a force that is directly proportional to the product of their masses and inversely proportional to the square of the distance between their centers.
### Mathematical Formulation
The gravitational force $\mathbf{F}$ between two masses $m_1$ and $m_2$ separated by a distance $r$ is given by:
$$
\mathbf{F} = -G \frac{m_1 m_2}{r^2} \hat{r}
$$
Where:
- $G$ is the gravitational constant ($6.674 \times 10^{-11} \ \mathrm{Nm^2/kg^2}$),
- $\hat{r}$ is the unit vector pointing from one mass to the other,
- $m_1$ and $m_2$ are the masses of the two bodies,
- $r$ is the distance between the centers of the masses.
### Gravitational Field
The gravitational field $\mathbf{g}$ at a distance $r$ from a mass $M$ is defined as the force per unit mass exerted by $M$:
$$
\mathbf{g} = -G \frac{M}{r^2} \hat{r}
$$
This field points toward the mass $M$ and has a magnitude:
$$
|\mathbf{g}| = G \frac{M}{r^2}
$$
#### Visualization: Gravitational Field
We visualize the gravitational field of a central mass $M$ in 2D slice.
```{python}
import numpy as np
import matplotlib.pyplot as plt
# Parameters
M = 2e22 # Mass of central body (kg)
G = 6.674e-11 # Gravitational constant
# Grid for visualization
x = np.linspace(-100, 100, 25)
y = np.linspace(-100, 100, 25)
X, Y = np.meshgrid(x, y)
R = np.sqrt(X**2 + Y**2)
# Avoid division by zero
R[R == 0] = np.nan
R[R <35] = np.nan
# Gravitational field
Fx = -G * M * X / R**3
Fy = -G * M * Y / R**3
# Plot field
plt.figure(figsize=(8, 8))
plt.quiver(X, Y, Fx, Fy, scale=1e10, color='blue')
plt.xlabel('x (m)')
plt.ylabel('y (m)')
plt.title('Gravitational Field')
plt.grid()
plt.gca().set_aspect('equal', adjustable='box')
plt.show()
```
### Orbital Motion
Gravitational force is responsible for the orbital motion of planets, moons, and satellites. The centripetal force required for circular motion is provided by gravity:
$$
F = \frac{mv^2}{r} = G \frac{mM}{r^2}
$$
From this, the orbital velocity $v$ of a body of mass $m$ around a central mass $M$ is:
$$
v = \sqrt{G \frac{M}{r}}
$$
The orbital period $T$ of the body can also be derived:
$$
T = 2\pi \sqrt{\frac{r^3}{GM}}
$$
<iframe src="html_anim/mechanics/gravity_orbits_1.html" width="100%" height="700px" frameborder="0"></iframe>
[html_anim/mechanics/gravity_orbits_1](html_anim/mechanics/gravity_orbits_1.html)
### Escape Velocity
Escape velocity is the minimum velocity required for an object to escape the gravitational pull of a planet or celestial body without further propulsion. It is given by:
$$
v_\text{escape} = \sqrt{2G \frac{M}{R}}
$$
where $M$ is the mass of the celestial body and $R$ is its radius.
### Numerical Example: Escape Velocity for Earth
- Mass of Earth ($M_\text{Earth}$): $5.972 \times 10^{24} \ \mathrm{kg}$
- Radius of Earth ($R_\text{Earth}$): $6.371 \times 10^6 \ \mathrm{m}$
```{python}
# Constants
G = 6.674e-11 # Gravitational constant (N m^2/kg^2)
M_earth = 5.972e24 # Mass of Earth (kg)
R_earth = 6.371e6 # Radius of Earth (m)
# Escape velocity calculation
v_escape = np.sqrt(2 * G * M_earth / R_earth)
print(f"Escape velocity for Earth: {v_escape/1000:.1f} km/s")
```
<iframe src="html_anim/mechanics/gravity_newton_canon.html" width="100%" height="700px" frameborder="0"></iframe>
[html_anim/mechanics/gravity_newton_canon](html_anim/mechanics/gravity_newton_canon.html)
### Gravitational Potential Energy
The gravitational potential energy $U$ of a two-body system is given by:
$$
U = -G \frac{m_1 m_2}{r}
$$
This energy is negative because the gravitational force is attractive, and the potential energy is zero at infinite separation.
### Applications of Universal Gravitation
1. **Planetary Orbits:** Predicting the motion of planets and moons in their orbits.
2. **Space Travel:** Calculating escape velocities and transfer orbits for spacecraft.
3. **Tidal Effects:** Understanding the gravitational interaction between Earth and Moon causing ocean tides.
4. **Astrophysics:** Studying gravitational interactions in galaxies and black holes.
### Small body around the Earth
```{python}
import numpy as np
import matplotlib.pyplot as plt
# Parametry
G = 6.674e-11 # Stała grawitacyjna (m^3/kg/s^2)
M = 5.972e24 # Masa większego ciała, np. Ziemi (kg)
m = 1000 # Masa mniejszego ciała, np. satelity (kg)
# Warunki początkowe
r0 = 7e6 # Początkowa odległość od środka masy większego ciała (m)
v0 = 5200 # Początkowa prędkość orbitalna (m/s)
theta0 = 0 # Początkowy kąt w radianach
# Czas
T = 6000 # Czas symulacji (s)
N = 100_000 # Liczba kroków czasowych
dt = T / N
# Tablice dla pozycji i prędkości
x = np.zeros(N)
y = np.zeros(N)
vx = np.zeros(N)
vy = np.zeros(N)
# Warunki początkowe
x[0] = r0 * np.cos(theta0)
y[0] = r0 * np.sin(theta0)
vx[0] = -v0 * np.sin(theta0)
vy[0] = v0 * np.cos(theta0)
# Numeryczna integracja (metoda Eulera)
for i in range(1, N):
r = np.sqrt(x[i - 1]**2 + y[i - 1]**2) # Odległość od środka masy
ax = -G * M * x[i - 1] / r**3 # Przyspieszenie w osi x
ay = -G * M * y[i - 1] / r**3 # Przyspieszenie w osi y
# Aktualizacja prędkości
vx[i] = vx[i - 1] + ax * dt
vy[i] = vy[i - 1] + ay * dt
# Aktualizacja pozycji
x[i] = x[i - 1] + vx[i - 1] * dt
y[i] = y[i - 1] + vy[i - 1] * dt
# Wizualizacja trajektorii
plt.figure(figsize=(8, 8))
plt.plot(x, y, label="Trajektoria")
plt.plot(0, 0, 'ro', label="Centralne ciało (np. Ziemia)")
plt.xlabel("x (m)")
plt.ylabel("y (m)")
plt.title("Trajektoria orbitalna małego ciała")
plt.legend(loc='lower right')
plt.grid()
plt.axis('equal')
plt.show()
```
### Earth-Moon System with small body
```{python}
import numpy as np
import matplotlib.pyplot as plt
# Parameters
G = 6.674e-11 # Gravitational constant (m^3/kg/s^2)
M_earth = 5.972e24 # Mass of Earth (kg)
M_moon = 7.348e22 # Mass of Moon (kg)
R_earth_moon = 3.844e8 # Distance between Earth and Moon (m)
m_probe = 1000 # Mass of the probe (kg)
# Initial positions
x_earth, y_earth = 0, 0
x_moon, y_moon = R_earth_moon, 0
# Launch parameters
launch_angle = 5 # Launch angle in degrees (relative to +x-axis)
initial_speed = 1500 # Initial speed of the probe (m/s)
# Convert launch angle to radians
launch_angle_rad = np.radians(launch_angle)
# Probe initial position and velocity
x_probe, y_probe = R_earth_moon / 2, 0 # Start halfway between Earth and Moon
vx_probe = initial_speed * np.cos(launch_angle_rad)
vy_probe = initial_speed * np.sin(launch_angle_rad)
# Time settings
T =15 * 24 * 3600 # Total simulation time (s)
N = 500000 # Number of time steps
dt = T / N # Time step (s)
# Arrays to store positions
x_positions = []
y_positions = []
# Simulation loop
for i in range(N):
# Distance between probe and Earth
r_earth = np.sqrt((x_probe - x_earth)**2 + (y_probe - y_earth)**2)
# Distance between probe and Moon
r_moon = np.sqrt((x_probe - x_moon)**2 + (y_probe - y_moon)**2)
# Gravitational accelerations
ax_earth = -G * M_earth * (x_probe - x_earth) / r_earth**3
ay_earth = -G * M_earth * (y_probe - y_earth) / r_earth**3
ax_moon = -G * M_moon * (x_probe - x_moon) / r_moon**3
ay_moon = -G * M_moon * (y_probe - y_moon) / r_moon**3
# Total acceleration
ax = ax_earth + ax_moon
ay = ay_earth + ay_moon
# Update velocities
vx_probe += ax * dt
vy_probe += ay * dt
# Update positions
x_probe += vx_probe * dt
y_probe += vy_probe * dt
# Store positions
x_positions.append(x_probe)
y_positions.append(y_probe)
# Convert to arrays
x_positions = np.array(x_positions)
y_positions = np.array(y_positions)
# Plot the results
plt.figure(figsize=(8, 8))
plt.plot(x_positions, y_positions, label="Probe Trajectory")
plt.plot(x_earth, y_earth, 'bo', label="Earth", markersize=10)
plt.plot(x_moon, y_moon, 'go', label="Moon", markersize=7)
plt.xlabel("x (m)")
plt.ylabel("y (m)")
plt.title("Trajectory of the Probe in the Earth-Moon System")
plt.legend()
plt.grid()
plt.axis('equal')
plt.show()
```
## Work and Kinetic Energy Relation
### Work
In three-dimensional motion, the work done by a force $\mathbf{F}$ acting on an object as it moves along a displacement $d\mathbf{r}$ is given by the **dot product**:
$$
dW = \mathbf{F} \cdot d\mathbf{r}
$$
Expanding in component form:
$$
dW = F_x dx + F_y dy + F_z dz
$$
Using **Newton’s second law** $\mathbf{F} = m \mathbf{a}$:
$$
dW = (m a_x dx + m a_y dy + m a_z dz)
$$
Since acceleration is the derivative of velocity:
$$
a_x = \frac{dv_x}{dt}, \quad a_y = \frac{dv_y}{dt}, \quad a_z = \frac{dv_z}{dt}
$$
And velocity components are:
$$
v_x = \frac{dx}{dt}, \quad v_y = \frac{dy}{dt}, \quad v_z = \frac{dz}{dt}
$$
Rewriting work in terms of velocity:
$$
dW = m v_x dv_x + m v_y dv_y + m v_z dv_z
$$
### Kinetic Energy
To determine the total work done from an initial velocity $\mathbf{v}_1 = (v_{1x}, v_{1y}, v_{1z})$ to a final velocity $\mathbf{v}_2 = (v_{2x}, v_{2y}, v_{2z})$:
$$
W = \int_{v_{1x}}^{v_{2x}} m v_x dv_x + \int_{v_{1y}}^{v_{2y}} m v_y dv_y + \int_{v_{1z}}^{v_{2z}} m v_z dv_z
$$
Solving each integral:
$$
W = m \left( \frac{v_{2x}^2}{2} - \frac{v_{1x}^2}{2} \right) + m \left( \frac{v_{2y}^2}{2} - \frac{v_{1y}^2}{2} \right) + m \left( \frac{v_{2z}^2}{2} - \frac{v_{1z}^2}{2} \right)
$$
Rewriting:
$$
W = \frac{1}{2} m \left( v_{2x}^2 + v_{2y}^2 + v_{2z}^2 \right) - \frac{1}{2} m \left( v_{1x}^2 + v_{1y}^2 + v_{1z}^2 \right)
$$
Since kinetic energy in 3D is defined as:
$$
KE = \frac{1}{2} m (v_x^2 + v_y^2 + v_z^2)
$$
we obtain:
$$
W = KE_2 - KE_1 = \Delta KE
$$
### Energy Conservation
In many physical situations, forces can be described by a **potential energy function**. A force $\mathbf{F}$ is called **conservative** if the work it does on an object moving from point $A$ to point $B$ is **independent of the path taken** and depends only on the initial and final positions. This allows us to define a scalar function, the **potential energy** $U$, such that:
$$
\mathbf{F} = - \nabla U
$$
This relation means that the force is the negative gradient of the potential energy function.
#### Examples of Potential Forces
There are several important cases of conservative forces where the **work-energy theorem** and **energy conservation** provide quick solutions to problems:
1. **Constant Force (Projectile Motion)**:
In uniform gravitational fields (e.g., near Earth's surface), the force is constant:
$$
\mathbf{F} = -mg \hat{j}
$$
The associated potential energy is given by:
$$
U = mg y
$$
This simplifies projectile motion calculations using energy conservation.
2. **Harmonic Oscillator (Spring Force)**:
A restoring force proportional to displacement, such as Hooke’s Law:
$$
\mathbf{F} = - k x \hat{i}
$$
The corresponding potential energy function is:
$$
U = \frac{1}{2} k x^2
$$
This is essential in analyzing oscillatory motion.
3. **Gravitational Force (Newtonian Gravity)**:
In a central gravitational field, the force obeys an inverse-square law:
$$
\mathbf{F} = - \frac{G m M}{r^2} \hat{r}
$$
The potential energy function is:
$$
U = -\frac{G m M}{r}
$$
This form is crucial for planetary motion and orbital mechanics.
### Energy Conservation
Since conservative forces allow potential energy to be defined, the **total mechanical energy** is conserved:
$$
E = KE + U = \text{constant}
$$
This principle allows for elegant solutions to motion problems without explicitly solving Newton’s second law.
#### Gravity example
$$
\frac{m v_1^2}{2}-\frac{G m M}{r_1}
=
\frac{m v_2^2}{2}-\frac{G m M}{r_2}
$$